Computationally Efficient Continuous-Time Model Predictive Control of a 2-DOF Helicopter via B-Spline Parameterization

This paper investigates one way to reduce the computational burden of continuous-time model predictive control (MPC) laws by representing the input/output signals and related models using B-spline functions. Such an approximation allows to implement the resulting feedback control law more efficiently, requiring less online computational effort. As a result, the proposed controller formulates the control signals as continuous polynomial spline functions. All constraints assumed over the prediction horizon are then expressed as constraints acting on the B-splines control polygon vertices. The performance of the proposed theoretical framework has been demonstrated with several real-time experiments using the well-known 2-DOF laboratory helicopter setup. The aim of the presented experiments was to track given step-like reference trajectories for pitch and yaw angles under notable parameter uncertainties. In order to suppress the influence of uncertainties, the control algorithm is implemented in an adaptive mode, equipped with the recursive least squares (RLS) estimation of model parameters and with the adaptation of stabilizing terminal set and terminal cost calculations. Thanks to the presented framework, it is possible to significantly reduce the computational burden, measured by the number of decision variables and input constrains, indicating the potential of the proposed concept for real-time applications, even when using embedded control hardware.


Introduction
The notion of spline has become an important and well-established tool in the modern approximation theory. At present, splines are successfully applied in an increasing number of scientific and technical disciplines, as can be seen, e.g., in [1][2][3]. Splines are, roughly speaking, functions which serve the purpose of approximating or modeling other functions or discrete data. Originally, splines were defined by I.J. Schoenberg in 1946 as univariate piecewise polynomial functions. One of the fundamental contributions to the theory of polynomial splines represents the discovery of compactly supported basis functions, the so-called B-splines, by Curry and Schoenberg in 1966. B-splines are piecewise polynomial functions which form a basis for the space of polynomial splines. Their importance arises from the fact that their properties make them highly suitable for practical computations. Being polynomial, they can be evaluated very quickly; being piecewise polynomial, they are very flexible. The number and profile of the B-spline functions is clearly determined by the distribution of the so-called knot points and it is possible to change their behavior by changing the position of these knots. The ability of splines to achieve a very high quality of approximation in various technical applications inspired the authors to use them when looking for alternative ways to reduce the computational complexity of continuous-time MPC laws. mentation, the developed adaptive B-spline-based CMPC controllers were deployed by means of the MATLAB/Simulink Desktop Real-Time software environment.

Spline and B-Spline Functions
In this section, we provide a brief summary of the basic properties of B-spline functions, conversions between B-and pp-spline representations, as well as their essential shape properties. For a comprehensive theory of splines, we refer the interested reader to [23][24][25][26]. In terms of content, the following mathematical background is mainly recalled from [16].

Basic Definitions
Let us first state some formal definitions. Definition 1. Let (ξ 0 =)0 < ξ 1 < . . . < ξ q < T x (= ξ q+1 ) be the subdivision of a closed finite-time interval [0, T x ] by q distinct (time) points. A function s(t), defined on the interval [0, T x ], is called a spline function of the order r > 0 (degree r − 1) and the defect d ef if the following two conditions hold: where α i denotes a multiplicity (or defect d ef ) of the knot ξ i . Let P r,ξ,α denote the linear space of spline functions for α = (α 1 , . . . α q ). If all interior knots, ξ i , are simple, i.e., α i = 1, i = 1, . . . , q, then P r,ξ,α = C r−2 [0, T x ]. In order to perform computations with splines, one must first choose a suitable representation, in which any member of P r,ξ,α can be written as a unique linear combination of properly chosen z basis functions such that Definition 1 is satisfied. A common choice is to use B-spline functions.

Definition 2.
A B-spline function N i,r,ξ (t) of order r > 0, with knots ξ i , . . . , ξ i+r , can be defined using the following recurrence relation: where the two fraction terms are interpreted as zero whenever ξ i+k − ξ i+1 = 0 and ξ i+k−1 − ξ i = 0, respectively.
From Definition 2, one can observe that N i,r,ξ (t), i = . . . 0, 1 . . . (i) have a local support, (ii) are positive on their supports and (iii) form a partition of unity. Every function s(t) satisfying Definition 1 then has a unique representation: with n(t) = [N 1,r,ξ (t), . . . , N z,r,ξ (t)] T , c = [c 1 , . . . , c z ] T , where N i,r,ξ (t) or shortly N i (t), i = 1, . . . , z, denote base functions of the spline space P r,ξ,α , and c i denotes the i-th B-spline coefficient of s(t). They are commonly referred to as control coefficients or control points, and the collection {c i } z i=1 of all control points is referred to as the control polygon of the spline. In the following, we will assume the spline represented as a linear combination of basis B-spline functions (2) as the approximation function. The spline design parameters thus are the following: • The order r and defect d ef of the spline; • The number q and location ξ of its knots; • The control coefficients (control polygon) {c i } z i=1 .

Conversion from B-to pp-Representation
A polynomial spline can be, by definition, written as where p i (t) are polynomial pieces or segments which represent the spline s(t) on interval [0, T x ]. The relation (3) is called a piecewise polynomial (pp-)representation of the spline s(t). Clearly, the pp-representation of spline s(t) is completely determined by the (q + 1)rdimensional vector of the polynomial coefficients p = [p 01 , . . . , p 0r , · · · , p q1 , . . . , p qr ] T . Given the B-representation (2) of the spline s(t), the vector p of its pp-representation can be easily computed according to p = P c, where rows of matrix P can be obtained by the differentiation of vector n T (t) in knots [23,27] for details. Given the vector p, the conversion from pp-representation to B-representation can be performed as follows: c = P −1 L p, which is more difficult because of the left inverse of matrix P; however, if it is a priori known that the approximated function lies in P r,ξ,α for a certain r, ξ, α, then P −1 L can be determined uniquely.

Shape Properties of B-Spline Curves
The fundamental shape properties of B-spline curves that we rely on in this article may be summarized using the following theorem. The reader is referred to, e.g., [26] for more details and proofs.
Theorem 1 (Shape properties of B-spline curves). Let s(t) be a B-spline curve of order r over the knot sequence ξ. Then, the following properties hold: In general, there is no endpoint interpolation; (ii) For ξ i ≤ t < ξ i+1 , s(t) lies in the convex hull of the r control points c i−r+1 , . . . , c i ; (iii) Local control: for t ∈ [ξ i , ξ i + 1], the curve is independent of c j for j < i − r + 1 and j > i; (iv) If r − 1 control points coincide, then the spline curve passes through this point and is tangent to the control polygon; (v) If r − 1 control points are on a line, then the spline curve touches this line; (vi) If r control points are on a line L, then s(t) ∈ L for ξ i ≤ t < ξ i+1 , i.e., an entire segment of the curve s(t) coincides with L; (vii) Its derivative is . . = ξ i+r−1 coincide, then s(t) = c i , i.e., the spline curve passes through a control point and is tangent to the control polygon.

Problem Formulation
The problem formulation will be presented using the example of the aforementioned laboratory 2-DOF helicopter, the hardware of which will be described later in Section 4. In some passages, it will follow formulations presented by the authors in [17] which led to the explicit solution to the spline-based CMPC problem.
Let us in general consider the laboratory 2-DOF helicopter as a continuous-time MIMO system described near a given operating point by the following linear time-invariant statespace form (in reality, in our application, it will be a two-input two-output (TITO) system): where A h , B h and C h are state-space matrices of appropriate dimensions and the elements of the state vector x h (t) ∈ R r u (r u denoted the order of the spline input signal) can be directly calculated from the spline derivatives of input/output signals u(t) ∈ R n u , y(t) ∈ R n y by using the multi-input single-output (MISO) technique formulated below.

Phase-Variable State-Space Description
Let us consider the MIMO system description (5) in the deterministic framework and apply it over the spline output and input signals, s y i (t) and s u j (t), i = 1 . . . n y , j = 1 . . . n u . Then, for the i-th spline output signal, s y i (t), of the system, the following MISO differential equation can be written: where coefficients α i,k and β l i,j , k = 1 . . . ρ α , l = 1 . . . ρ β , are real constant scalars. For the degrees ρ α and ρ β , it holds that ρ α ≤ (r y − 1) and ρ β ≤ (r u − 1), respectively. Having specified spline input-output derivatives, by applying a common technique of moving data windows it is easy to estimate the vector of unknown parameters of (6) through linear regression: and a linear square approach where i s (t) denotes an equation error. Several methods are available for finding possible state-space forms from the inputoutput map (6). The well-known phase-variable canonical forms, see, e.g., [28], are particularly useful for our purposes because state variables can be defined in terms of the plant variables and their derivatives. With this, we can choose the ρ α -dimensional state T , of a phase-variable canonical form related to the i-th spline output, as with (ρ α × ρ α ) square matrices S y i , S u i,j defined as The foregoing equations can be arranged in the following MISO state-space form: where By assuming i = 1 . . . n y , the set of MISO realizations (9) can be used to create the following simple MIMO state-space realization of the system which is valid near the given operating point:ẋ with order n = n y .ρ α , and the state and output vectors defined as It is easy to verify that The realization (10) is a straightforward connection of canonical realizations (9) for each system output and therefore is not the minimal one. Nevertheless, in this realization, the n-dimensional state vector x m (t) is accessible to direct measurement.
If the system's operating point varies in time-and this is our case-then the description (10) also becomes time-variant and can be written as: where the matrices A m (t) and B m (t) can be easily identified online using, e.g., algorithms of the recursive least squares method; thus, the overall design of the control scheme will enter the adaptive mode.

B-Spline Parameterization of CMPC Formulation
In this subsection, we discuss the CMPC formulation from the viewpoint of the application of B-spline parameterization as a starting point for obtaining computationally efficient control algorithms. In order to be able to introduce the integral action in the resulting formulations, we are further more interested in the derivatives of control input signals contained in the vectoru(τ). Using the expansion (2) with the same choice of Bspline function parameters for all inputs (i.e., u i (τ) ≡ s u i (τ) = n(τ) T c u i , u i (τ) ∈ P r u ,ξ,α , i = 1, . . . , n u ), these can be obtained asu . . , n u andṅ(τ) denoting a vector of the first derivatives of the used B-spline functions.
By taking the derivative of both sides of (11) with u(t) ≡ s u (t) and y(t) ≈ s y (t), the state-space model can be rewritten in the following augmented form: where Assuming a time window given by τ ∈ [t k , t k + T h ], where t k is the current time and T h denotes the prediction horizon, the predicted state at time t k + τ can be obtained by solving the differential Equation (13a) as follows: where for simplicity we used the notation A k = A(t k ) and B k = B(t k ). Now, if we employ the B-spline functions-based expansion (12) to approximate the control input derivatives, we can substitute them into the prediction Equation (14), which thus becomes parameterized in c u : The objective of the velocity-form CMPC design considered in this work is to drive the predicted system outputs, y(τ), as close as possible to predefined reference trajectories y s (τ), ideally to solve the following infinite-horizon control problem: where within the quadratic objective (16a), Q y (τ) 0 and Q u (τ) 0 denote the output and input weighting matrix, respectively. In the constraints (16b) and (16c), (u min , u max ) and (u min ,u max ) denote bounds imposed on control signals and their derivatives, respectively. Now, applying the quasi-infinite horizon approach to guarantee the closed-loop stability, our goal is to solve the following finite-horizon continuous-time MPC problem: where Q x (τ) = C T Q y (τ)C, x s (τ) denotes the state reference trajectories defined by the setpoints, andx(t k + T h |t k ) denotes a deviation from the steady-state target calculated for the current reference state x s (τ) at τ = T h . In addition, the stability and recursive feasibility of the control problem are ensured by assuming a terminal cost weighted with Q h,k 0 in (17a), and by assuming a terminal set Ω k constraint (17d). This set has to be invariant under a local linear state feedback u = Fx, virtually acting for τ ∈ [t k + T h , ∞[ and feasible with (17b) and (17c). For the calculation of Ω k and Q h,k , a simple computational procedure proposed in [29] can be adopted. The adaptive CMPC problem (17) can be tackled using B-spline parameterization which leads to the following formulation: The reformulated objective (18a) was obtained by substituting (12) and (15) into (17a), while expressing the setpoints as spline functions using the relation (2) as follows: . . , n y , and c s (t k ) is a vector of the splines' control points expressing the trajectories of setpoints prescribed by the user for the desired system outputs' behavior over the prediction horizon τ ∈ [t k , t k + T h ]. The matrices H k , G k and Q h,k can be calculated, using numerical integration, as follows: Note that the continuous-time constraints (17b) and (17c) had to be reformulated to suitable finite-dimensional forms given by (18b) and (18c), respectively. These must guarantee the intersample behavior of the spline control signal, which is achieved by a proper bounding of its control polygon c u (k). Taking into account the basic shape properties of B-splines, (18b) represents the amplitude constraints (17b) of the control signal u(t) imposed over a horizon [t k , t k + T h ], where c min u and c max u are the min and max values of the spline control coefficients c u (k) computed using the B-spline approximation of the given u min and u max values. Similarly, (18c) approximates the constraints (17c) on the derivative of the control signalu(t), where c min u,∆ and c max u,∆ are the min and max values of the differences in the spline control coefficients, which can, together with entries of matrix A ∆ , be computed based on the location of the spline knots. Next, (18d) represents the terminal constraints transformed from (17d) using the prediction model (15), where denote the vectors bounding the terminal set Ω k . In order to keep the spline function s u i (t) in the space P r u ,ξ,α for the given (r u , ξ, α), the equality constraints (18e) and (18f) are used to enforce continuity between the implemented and the projected spline control signals at the beginning of the prediction horizon [t k , t k + T h ]. Finally, the constraint (18g) is added to improve the stability of the control signals by requiring a zero derivative of the projected spline control signal at the end of the prediction horizon [t k , t k + T h ].
From the optimization perspective, the adaptive spline-based stabilizing CMPC problem in velocity form, (18), is a quadratic program (QP), which can be recast in the following simplified form: Solving QP (19) implicitly for a current state x(t k ) and given setpoints c s (t k ) yields a vector of optimal control coefficients c u . . , n u . According to [15,16], c u (k) is subsequently converted to piecewise polynomial (pp)-representation of the individual optimal spline control signals s u i (τ), τ ∈ [t k , t k + T h ], i = 1, . . . , n u given by their polynomial coefficients contained in a vector p (k) = P c u (k).
The optimization problem (18) implies that we are looking at the control signals u i (t) from the viewpoint of a selected distance T between the knot points of the spline function s u i (t). In real-time implementation, the distance T corresponds to the control period, in which the first polynomial segments of the respective optimal spline input signals, i.e., , are applied for control. This is usually performed with a much shorter implementation period T g T; see Figure 1 for illustration. This means that we calculate the values of the control signal s u i (t k + τ) for the relative time variable τ = jT g , j = 1, . . . , n g , T = n g T g , i = 1, . . . , n u and implement them for control using a common zero-order hold. The parameter n g is user-defined and its choice is important because it determines the degree of continuity of the generated signal. The optimization problem (18) implies that we are looking at the control signals u i (t) 209 from the viewpoint of a selected distance T between the knot points of the spline function 210 s u i (t). In real-time implementation, the distance T corresponds to the control period, in 211 which the first polynomial segments of respective optimal spline input signals, i.e. s ⋆ u i (τ), 212 τ ∈ [t k , t k + T], are applied for control. This is usually performed with a much shorter-213 implementation period T g ≪ T; see Figure 1 for illustration. This means that we calculate 214 the values of control signal s ⋆ u i (t k + τ) for the relative time variable τ = jT g , j = 1, . . . , n g , 215 T = n g T g , i = 1, . . . , n u and implement them for control using a common zero-order hold. 216 The parameter n g is user-defined and its choice is important because it determines the 217 degree of continuity of the generated signal. Note that during the prediction horizon [t k , t k + T h ] the distance T between the knots 220 of the projected spline input signal s u i (τ), τ ∈ [t k , t k + T h ] can be chosen as βT, β = 1, 2, . . ., 221 introducing a property comparable with the well-known move blocking techniques used 222 in MPC.

223
This is an important design parameter which-in conjunction with spline continuity 224 conditions and the spline defect-allows us to tune various controller setups. Essentially, 225 the number q of polynomial segments of the projected spline control signal determines the 226 degrees of freedom for a given task. A higher value of q provides a more active control 227 signal while a lower value leads to control which is more smooth and sluggish. Note that 228 small values of q, q ≥ 1 significantly reduce the QP computational effort.

229
The projected polynomial segments are interconnected in interior knot points of the 230 sequence ξ h which spans over the prediction horizon, [t k , t k + T h ]. In general, the place-231 ment of these knots may be considered as another design tool. As mentioned earlier, we 232 prefer a uniform location of these knots with a selected distance T u , taken as T u = βT. 233 The distance then becomes a remarkable tuning knob for setting the overall length T h of 234 the prediction horizon under the same number q of its segments. Using the distance T u 235 as a variable, we can extend or shorten the prediction horizon satisfying the same compu-236 tational burden due to the same dimension of the decision vector c u . The chosen length 237 of T u considerably dominates activity of the projected spline control signal. Consequently, 238 there is no need for any control weighting in the cost and it is possible to set w u (t) = 0, 239 even for control of non-minimum phase plants. The setting w u (t) = 0 was assumed in all 240 experiments reported in this work. 241 Figure 1. Illustration of how a spline control signal is generated. • denotes a knot point, T-sampling period, T g -implementation (generation) period.

Number of Polynomial Segments of Projected Spline Control Signals
Note that during the prediction horizon [t k , t k + T h ] the distance T between the knots of the projected spline input signal s u i (τ), τ ∈ [t k , t k + T h ] can be chosen as βT, β = 1, 2, . . ., introducing a property comparable with the well-known move-blocking techniques used in MPC. This is an important design parameter which-in conjunction with spline continuity conditions and the spline defect-allows us to tune various controller setups. Essentially, the number q of the polynomial segments in the projected spline control signal determines the degrees of freedom for a given task. A higher value of q provides a more active control signal, while a lower value leads to control which is more smooth and sluggish. Note that small values of q, q ≥ 1 significantly reduce the QP computational effort.
The projected polynomial segments are interconnected in the interior knot points of the sequence ξ h which spans over the prediction horizon, [t k , t k + T h ]. In general, the placement of these knots may be considered as another design tool. As mentioned earlier, we prefer a uniform location of these knots with a selected distance T u , taken as T u = βT. The distance then becomes a remarkable tuning knob for setting the overall length T h of the prediction horizon under the same number q of its segments. Using the distance T u as a variable, we can extend or shorten the prediction horizon satisfying the same computational burden due to the same dimension of the decision vector c u . The chosen length of T u considerably dominates the activity of the projected spline control signal. Consequently, there is no need for any control weighting in the cost and it is possible to set w u (t) = 0, even for control of non-minimum phase plants. The setting w u (t) = 0 was assumed in all the experiments reported in this work.
One may tune different controllers via alternating continuity conditions and the spline defect in selected knots of the projected spline control signal. For example, it is possible to make-like in generalized predictive control-an assumption of zero increments of the projected control signal, starting from a certain time instant in the horizon, and thereby to impose an actual control horizon within the given prediction horizon. This can be simply achieved by constraining the knot derivatives of a certain number of last polynomial segments to zero while setting spline defects in the involved joining knots to zero as well. Applying the above procedure to a chosen number of segments, we can set a proper length of the control horizon, see, e.g., Figure 2a,b. For illustrative purposes, in these figures, two arrangements of third-order polynomial spline segments within the prediction horizon are presented. In Figure 2a, the derivatives of the last segment are set to zero, while in a similar arrangement in Figure 2b, a zero defect in the next-to-last knot is added. The impact on the actual control horizon length is evident. One may tune different controllers via alternating continuity conditions and the spline 242 defect in selected knots of the projected spline control signal. For example, it is possible 243 to make-like in generalized predictive control-an assumption of zero increments of the 244 projected control signal, starting from a certain time instant in the horizon, and thereby 245 to impose an actual control horizon within the given prediction horizon. This can be sim-246 ply achieved by constraining the knot derivatives of a certain number of last polynomial 247 segments to zero while setting spline defects in the involved joining knots to zero as well. 248 Applying the above procedure to chosen number of segments we can set a proper length 249 of the control horizon, see e.g. Figure 2a and Figure 2b. For illustrative purposes in these 250 figures two arrangements of third-order polynomial spline segments within the prediction 251 horizon are presented. In Figure 2a derivatives of the last segment are set to zero, while 252 in a similar arrangement in Figure 2b a zero defect in the next-to-last knot is added. The 253 impact on the actual control horizon length is evident.

254
Following the receding horizon strategy, only that piece of the first polynomial seg-255 ment of the spline control signal which corresponds to the first sampling period T is really 256 applied to the controlled system as per Figure 1. For current plant data y k , c u (k), and reference values c s (t k ), using a recursive least 260 square algorithm update the model (10) (19) to obtain the vector of optimal 265 control coefficients c ⋆ u (k); which is then used to 266 5. calculate the vector p ⋆ (k) containing polynomial coefficients of all segments, using 267 the relation (4) as p ⋆ (k) = P c ⋆ u (k); and 268 6. following the receding horizon strategy apply with the possibly shortest implemen-269 tation period T g , T g ≪ T the first polynomial segments of p ⋆ (k) to control the sys-270 tem during the time interval [t k , t k + T], which means calculate the control signals 271 s ⋆ u i (t k + τ) for a relative time variable τ = jT g , j = 1, . . . , n g , T = n g T g , i = 1, . . . , n u , 272 as spline polynomial segments of order r u , and implement them for control using a 273 common zero-order hold; and finally 274 7. repeat the procedure from step 1. for the next sampling instant t k+1 = t k + T.
275 Figure 6 shows a simplified block scheme of the above algorithm. Following the receding horizon strategy, only that piece of the first polynomial segment of the spline control signal which corresponds to the first sampling period T is really applied to the controlled system, as per Figure 1.
Overall, the proposed adaptive B-spline-based stabilizing CMPC algorithm boils down to performing the following steps at each time instant t k : 1. For current plant data y k , c u (k), and reference values c s (t k ), using a recursive least squares algorithm update the model (10) and calculate the state vector x(t k ); use it to 2. adapt the calculation of bounds x min Ω k and x max Ω k of the terminal set Ω k , and terminal cost Q h,k ; and to 3. update the matrices H k , G k and Q s,k ; 4. solve the optimal control problem (18) as the QP (19) to obtain the vector of optimal control coefficients c u (k); which is then used to 5. calculate the vector p (k) containing polynomial coefficients of all segments, using the relation (4) as p (k) = P c u (k); and 6. following the receding horizon strategy apply with the possibly shortest implementation period T g , T g T the first polynomial segments of p (k) to control the system during the time interval [t k , t k + T], which means calculate the control signals s u i (t k + τ) for a relative time variable τ = jT g , j = 1, . . . , n g , T = n g T g , i = 1, . . . , n u , as spline polynomial segments of order r u , and implement them for control using a common zero-order hold; and finally 7. repeat the procedure from step 1. for the next sampling instant t k+1 = t k + T. Figure 3 shows a simplified block scheme of the above algorithm.    Our laboratory helicopter physical model [7] consists of a base and a beam carrying 287 at its ends two propellers-the main and the tail one-driven by DC motors; see Figure 4. 288 The beam is attached to its base via an articulated joint which allows the beam to rotate so 289 that its ends move on spherical surfaces. This connection enables two degrees of freedom 290 of the helicopter body movement-rotation around the horizontal axis, described by the 291 pitch angle, and rotation around the vertical axis, described by the yaw angle. Both angles 292 are measured using high-resolution incremental encoders. The model was designed at the 293 authors' workplace and represents a low-cost alternative to the commercially available 294 models, such as [30], except for its limited rotation around the vertical axis imposed by 295 power and data cables.

296
The axes of the main and the tail propeller as well as the vertical and the horizontal 297 helicopter axis are usually perpendicular to each other so that the movement in the verti-298 cal plane and the movement in horizontal plane are each affected by the thrust of only one 299 propeller. A counter-weight fixed to the beam determines a stable equilibrium position. 300 The system is balanced in such a way that when the motors are not powered, the main 301 propeller end of the beam is lowered. As it is usual for similar small-scale laboratory he-302 licopter setups, the control is achieved exclusively by controlling speeds of the propellers 303 at a fixed angle of attack. The range of the helicopter body rotation is ±30 • in pitch angle 304 and ±145 • in yaw angle. The plane of our main propeller is also slightly (∼2 • ) deviated 305 from the horizontal one to strengthen the coupling effect with the tail propeller. The main 306

Experimental Results and Discussion
This section describes the experiments carried out to verify the functionality and efficiency of the proposed concept of the B-spline CMPC algorithm. In the following subsections, we briefly explain the experimental 2-DOF helicopter platform and software tools, and after that, we illustrate some of the obtained results on four examples of tracking the predefined reference trajectories for the pitch and yaw angles. Similar platforms are well-known in the literature and are commonly used as benchmarks for testing various control techniques. In our case, the platform allows to easily demonstrate the basic features of the design and performance of the proposed B-spline-based CMPC controller.

Experimental Setup
Our laboratory helicopter physical model [7] consists of a base and a beam carrying at its ends two propellers-the main and the tail one-driven by DC motors; see Figure 4. The beam is attached to its base via an articulated joint which allows the beam to rotate so that its ends move on spherical surfaces. This connection enables two degrees of freedom of the helicopter body movement-rotation around the horizontal axis, described by the pitch angle, and rotation around the vertical axis, described by the yaw angle. Both angles are measured using high-resolution incremental encoders. The model was designed at the authors' workplace and represents a low-cost alternative to the commercially available models, such as [30], except for its limited rotation around the vertical axis imposed by power and data cables.
The axes of the main and the tail propellers as well as the vertical and the horizontal helicopter axes are usually perpendicular to each other so that the movement in the vertical plane and the movement in the horizontal plane are each affected by the thrust of only one propeller. A counter-weight fixed to the beam determines a stable equilibrium position. The system is balanced in such a way that when the motors are not powered, the main propeller end of the beam is lowered. As is usual for similar small-scale laboratory helicopter setups, the control is achieved exclusively by controlling the speeds of the propellers at a fixed angle of attack. The range of the helicopter body rotation is ±30°in the pitch angle and ±145°in the yaw angle. The plane of our main propeller is also slightly (∼2°) deviated from the horizontal one to strengthen the coupling effect with the tail propeller. The main source of random uncertainties in our laboratory helicopter design is the persistent vibration of the tail endpoint of the beam, which manifests itself in the increased fluctuation of the measured value of the helicopter yaw angle during the experiments.  source of random uncertainties in our laboratory helicopter design is the persistent vibra-307 tion of the tail end point of the beam, which manifests itself in the increased fluctuation of 308 the measured value of the helicopter yaw angle during experiments.

309
The described laboratory helicopter model can be represented as a nonlinear multi-310 variable system with two inputs:

317
The interface voltages u 1 and u 2 applied for setting of the helicopter inputs are con-318 verted to appropriate voltage values that drive the propeller motors. The output y 1 de-319 notes the pitch angle in the vertical plane between the longitudinal axis of the helicopter 320 body and the horizontal axis, and the output y 2 denotes the yaw angle in the horizontal 321 plane between the longitudinal axis of the helicopter body and its zero (initial) position.

322
The voltage driving the main motor and the voltage driving the tail motor affect both 323 the pitch angle and the yaw angle, hence we can say that the interaction makes the system 324 multivariable. It is also worth mentioning that the used mechanical simplification with 325 fixed-angle propeller blades does not necessarily translate into simplified dynamics. On 326 the contrary, the input torques and forces are applied via aerodynamical effects, as well 327 as additional coupling effects appearing between the helicopter body and propellers dy-328 namics, due to the reaction forces and torques arising at the acceleration or deceleration of 329 the propellers. These increased cross-coupling effects have important implications on the 330 control of the helicopter system and make its dynamics, as already mentioned, partially 331 uncertain and time-variable. This is the reason why we consider the design of control in 332 an adaptive mode with RLS estimation of the model parameters and adaptation of calcu-333 lation of the stabilizing terminal set. To solve the problem of RLS estimator wind-up, the 334 estimator with directional forgetting from [21] was implemented. This estimator ensures 335 the convergence of the estimation and avoids large changes in model parameters. The described laboratory helicopter model can be represented as a nonlinear multivariable system with two inputs: • u 1 -voltage driving propeller speed of the main motor; • u 2 -voltage driving propeller speed of the tail motor.
These are manipulated in the interface range of 0 V to 10 V. There are also two outputs: • y 1 -pitch (elevation) angle; • y 2 -yaw (azimuth) angle.
These are measured in degrees. The interface voltages u 1 and u 2 applied for setting the helicopter inputs are converted to appropriate voltage values that drive the propeller motors. The output y 1 denotes the pitch angle in the vertical plane between the longitudinal axis of the helicopter body and the horizontal axis, and the output y 2 denotes the yaw angle in the horizontal plane between the longitudinal axis of the helicopter body and its zero (initial) position.
The voltage driving the main motor and the voltage driving the tail motor affect both the pitch angle and the yaw angle; hence, we can say that the interaction makes the system multivariable. It is also worth mentioning that the used mechanical simplification with fixed-angle propeller blades does not necessarily translate into simplified dynamics. On the contrary, the input torques and forces are applied via aerodynamical effects, as well as additional coupling effects appearing between the helicopter body and propellers dynamics, due to the reaction forces and torques arising at the acceleration or deceleration of the propellers. These increased cross-coupling effects have important implications on the control of the helicopter system and make its dynamics, as already mentioned, partially uncertain and time-variable. This is the reason why we consider the design of the control in an adaptive mode with RLS estimation of the model parameters and the adaptation of the calculation of the stabilizing terminal set. To solve the problem of RLS estimator wind-up, the estimator with directional forgetting from [21] was implemented. This estimator ensures the convergence of the estimation and avoids large changes in the model parameters.
The proposed spline-based CMPC scheme to control the laboratory helicopter shown in Figure 4 was implemented using the MATLAB / Simulink Desktop Real-Time prototyping suite on a mini PC equipped with a 2.8 GHz CPU and 16 GB of RAM. Communication with the helicopter testbed was handled using the Humusoft's MF644 multifunction desk-top I/O card. The underlying QP (19) of the MPC problem is solved repeatedly at each control instant T using a parametric active-set algorithm implemented in the open-source software package qpOASES [31].

Experimental Results
In order to demonstrate the functionality and computational efficiency of the proposed spline-based CMPC scheme, we present four experiments assuming values of the design parameters listed in Table 1. These parameters are complemented by the selection of weighing matrices in the objective function (18), which were set as follows: • As outlined, the input weighing matrix Q u was set as zero for all experiments. • The output weighing matrix Q y was chosen as an exponential type with different exponential factors, λ 1 = 0.995 for the pitch angle and λ 2 = 0.967 for the yaw angle. The parameters of the adaptive directional forgetting used in the experiments were chosen as follows: • The minimum value of the adaptive forgetting factor was set to 0.97. • The expected value of the adaptive forgetting factor in a steady state was set to 0.99.  All four controllers were tested to follow the same sequence of nonsimultaneous step changes in the pitch angle and yaw angle references. The experiments were conducted for a selection of two different implementation periods, namely T g = 0.04 s (experiments E1 and E2) and T g = 0.02 s (experiments E3 and E4). In doing so, the chosen setting of the quality parameters of the B-spline representation (spline orders, number and location of interior knots) in these experiments led to the number n u (r u + q) of decision variables (i.e., the length of the optimizer vector) varying in the range of 10 to 14, as listed in Table 1.
The obtained tracking results for experiments E1 and E2 are presented in Figure 5a. One may notice the quite expected consequence that increasing the length of the optimizer improves the overall tracking performance, which is also evidenced by means of standard deviations of the tracking errors listed in Table 2. The comparison of the experiments is also supplemented by the profiles of the applied spline input signals, see Figure 5b, respecting the constraints imposed on the amplitude and the derivative. Finally, Figure 5c shows the numbers of iterations of the employed QP solver. Its addition is merely illustrative but to an extent confirms that lengthening the optimizer leads (under the same conditions) to lower numbers of required iterations.
The experiments were repeated with the implementation period T g = 0.02 s (experiments E3 and E4). As evidenced by Figure 6a and Table 1, the shorter implementation period contributed to an increase in the quality of tracking. Figure 6b,c illustrate the corresponding profiles of the spline input signals and solver iterations, respectively. We can also observe from the performed experiments that the quality of tracking the reference trajectories was significantly affected by the existing cross-couplings, characteristic for the controlled helicopter system. This is owed mainly to the mutual link between the pitch angle and the yaw angle, which in this state of the solution we tried to suppress by the mentioned selection of the output weighing matrix Q y .
We also remark that the zero input weighing matrix sufficed to ensure a stabilizing control. This observation is a by-product of the conditions (18e)-(18g) imposed on the continuity and zero derivative of the splines at the end of the prediction horizon.
In terms of evaluating the computational complexity, the main intention of this article is not to compare the proposed approach with other MPC formulations. Nevertheless, if we want to assess the computational burden of the proposed B-spline-based CMPC procedure compared with the standard discrete-time MPC procedure virtually acting on the same prediction horizon, we can approach it by comparing the sizes of the parameterizations that would be required by these procedures within the implementation, see Tables 3-5. They show that the B-spline parameterization is markedly leaner, i.e., computationally more efficient than the parameterization of a comparable standard discrete-time formulation, either in the case of parameterization assuming faster sampling given by control period T g ( Table 5) or slower sampling given by T (Table 4). Regarding these comparisons, we may also point out that the CMPC procedure in addition enables to guarantee the intersample behavior-a property that may be approximated by the discrete-time formulation with shorter sampling time T g , presented in Table 5.     300T g 600 2400 E3 150T g 300 1200 E4 200T g 400 1600

Conclusions
This paper presented a numerically efficient approach to the synthesis of an adaptive velocity integration version of a continuous-time MPC based on B-spline parameterization. The design procedure is shown for multivariate systems represented by a set of MISO models. The implementation itself was carried out in an adaptive mode, with the adaptation of the model and the stabilizing terminal set, in order to improve the tracking capabilities of the controlled system in an uncertain and stochastic environment.
The resulting control strategy was applied for the pitch and yaw angle control of a 2-DOF laboratory helicopter. The performed real-time experiments demonstrated that the proposed B-spline CMPC by design represents a computationally viable alternative to the standard discrete-time MPC implementation. The obtained results are documented both graphically and numerically by capturing the standard deviations of the control errors. At the same time, they allow the reader or user to clearly understand the impact of optional parameters of the B-spline parameterization on the resulting tracking performance.
Finally, we remark that the proposed algorithm is typically coded on a single processor, although a dual-processor implementation would enable to separate the generation of the input polynomial segment from the optimization process, which could further contribute to the reduction in the tracking error.